# Load necessary libraries #
library(devtools)
library(SlicerMorphR)
library(rgl)
library(geomorph)
library(readxl)
library(ks)
# Set your working directory
setwd("INSERT_PATH_HERE")
##############################################
# Load SlicerMorph Data and Preprocessing #
##############################################
# Choose the Log File from Slicer
# (Ensure a line break is added at the end of the .log file before loading)
# Choose the SlicerMorph .log file
SM.log.file <- file.choose()
SM.log <- parser(SM.log.file, forceLPS = TRUE)
# Extract Landmark Data
Data <- SM.log$LM
# Create geomorph dataframe
datagdf <- geomorph.data.frame(landmarks = Data)
# Perform Procrustes Superimposition
Pcoords <- gpagen(datagdf$landmarks)
# Save Procrustes Coordinates
save(Pcoords, file = "filename.bin")
##############################################
# Semi-Landmark Sliding #
##############################################
# Extract Semi-Landmarks Data
Data2 <- SM.log$semiLMs
# Sliding using **bending energy minimization**
PcoordsSlid <- gpagen(A = Data, surfaces = as.numeric(Data2), ProcD = FALSE, print.progress = TRUE)
save(PcoordsSlid, file = "filenameSlidLandmarks.bin")
# Extract and Save Centroid Size (Csize)
write.table(PcoordsSlid$Csize, file = "FileName1.tsv", sep = "\t")
##############################################
# Merge with Master Dataset #
##############################################
# Load Pcoords Lists
PcoordsList <- read.table(file = 'FileName1.tsv', sep = '\t', header = TRUE)
# Load Master Dataset
Masterfile <- read_excel("path/to/YourMasterfile.xlsx")
# Merge Pcoords with Master Dataset
Table_SSL <- merge(PcoordsList, Masterfile, by = "LabID", sort = FALSE)
# Save Merged Tables
write.csv(Table_SSL, file = "FileName2.csv", row.names = FALSE)
##############################################
# Create and Save Final Analysis Dataset #
##############################################
# SSL Dataset
spec <- read.csv("FileName2.csv", header = TRUE)
YourName.dt <- list(gpa.sh = PcoordsSlid$coords, cs = PcoordsSlid$Csize, spec = spec)
save(YourName.dt, file = "filename.dt.bin")
##############################################
# Symmetry Analysis #
##############################################
## Define Landmark Pairs for Symmetry Analysis - Replace by your own paired landmarks
lm.pairs <- matrix(c(1, 2, 4, 5, 8, 9, 11:20, 22:85, 87:128, 130:177, 179, 180, 182, 183, 185, 186, 188:301, 304:399, 401:412, 414:427, 429:446, 448:461, 465:600, 602:625, 627:652, 656:665, 667:694, 696:701), ncol = 2, byrow = T)
## Perform Symmetry Analysis
asym <- bilat.symmetry(PcoordsSlid$coords, ind = dimnames(PcoordsSlid$coords)[[3]],
object.sym = TRUE, land.pairs = lm.pairs, iter = 9)
symm.sh <- asym$symm.shape
# Save Symmetry Data
YourNewName.dt <- list(gpa.sh = PcoordsSlid$coords, cs = PcoordsSlid$Csize,
symm.sh = symm.sh, spec = spec)
save(YourNewName.dt, file = "YourNewfilename.dt.bin")
# Load necessary libraries
library(geomorph)
library(readr)
library(ks)
library(rgl)
# Set working directory (modify as needed)
setwd("INSERT_PATH_HERE")
# Load dataset
load("YourNewfilename.dt.bin")
######################################################################
# PCA Analysis #
######################################################################
# Perform Principal Component Analysis (PCA)
YourNamePCA <- gm.prcomp(YourNewfileName.dt$symm.sh)
# Access eigenvalues
eigenvalues <- YourNamePCA$sdev^2
# Calculate proportion of variance explained by each principal component
variance_explained <- eigenvalues / sum(eigenvalues)
# Print proportion of variance explained by each principal component
print(variance_explained)
######################################################################
# 3D PCA Plot #
######################################################################
# Plot - Replace x by your number of groups in the categorical variable
cols <- rainbow(x)
plot3d(YourNamePCA$x[,1:3], col = cols[as.factor(YourNewfileName.dt$spec$YourGroup)], size = 20)
# Add legend
legend3d("topright", legend = unique(as.factor(YourNewfileName.dt$spec$YourGroup)), col=c(unique(cols[as.factor(YourNewfileName.dt$spec$YourGroup)])), pch = 10)
# Add text
text3d(YourNamePCA$x[,1:3], texts = YourNewfileName.dt$spec$LabID, pos = 1, cex = 0.8)
######################################################################
# 3D Cloud Plot #
######################################################################
# Extract PCA scores
Score <- YourNamePCA$x[,1:3]
Labels <- factor(x=YourNamePCA$YourGroup)
# Kernel density estimation for group classification
Hscv1.Bd <- Hkda(x = YourNamePCA$x[,1:3], x.group = as.factor(YourNewfileName.dt$spec$YourGroup), bw = "scv", pre = "sphere")
# Perform kernel discriminant analysis
Tt.kda3.Bd <- kda(x = YourNamePCA$x[,1:3], x.group = as.factor(YourNewfileName.dt$spec$YourGroup), Hs = Hscv1.Bd)
# Plot 3D Cloud
plot(Tt.kda3.Bd, colors = c("YourColour","YourColour","YourColour", "etc"), drawpoints=TRUE, col.pt=c("YourColour","YourColour","YourColour", "etc"), box=FALSE, display="rgl")
# Add text
text3d(YourNamePCA$x[,1:3], texts = YourNewfileName.dt$spec$LabID, pos = 1, cex = 0.8)
# Perform PERMANOVA IN PAST
# Export the Pc Scores #
write.table(YourNamePCA$x, file="FileName3.tsv", sep="\t") # tsv because of space
######################################################################
# Plot 3D Coastal Only #
######################################################################
# Filter out group of specimens (e.g., filter out offshore)
# Identify indices for specimens we want to keep
Indice <- YourNewfileName.dt$spec$YourGroup != "offshore"
# Subset the dataset based on offshore specimens
YourNewfileName.dt$symm.sh <- YourNewfileName.dt$symm.sh[, , Indice, drop = FALSE]
YourNewfileName.dt$cs <- YourNewfileName.dt$cs[Indice, drop = FALSE]
YourNewfileName.dt$spec <- YourNewfileName.dt$spec[Indice, ]
# Create Geomorph Data Frame (GDF) without offshore specimens
gdfNoOffshore <- geomorph.data.frame(
shape = YourNewfileName.dt$symm.sh,
ecotype = YourNewfileName.dt$spec$YourGroup,
ind = YourNewfileName.dt$spec$LabID,
Csize = YourNewfileName.dt$cs
)
# Perform PCA
YourNamePCA <- gm.prcomp(YourNewfileName.dt$symm.sh)
######################################################################
# 3D PCA Plot #
######################################################################
# Replace x by your number of group in the categorical variable
cols <- rainbow(x)
plot3d(YourNamePCA$x[,1:3], col = cols[as.factor(YourNewfileName.dt$spec$YourGroup)], size = 20)
# Add a legend
legend3d("topright", legend = unique(as.factor(YourNewfileName.dt$spec$YourGroup)), col=c(unique(cols[as.factor(YourNewfileName.dt$spec$YourGroup)])), pch = 10)
# Add text
text3d(YourNamePCA$x[,1:3], texts = YourNewfileName.dt$spec$LabID, pos = 1, cex = 0.8)
######################################################################
# 3D Cloud Plot #
######################################################################
# Extract PCA scores
Score <- YourNamePCA$x[,1:3]
Labels <- factor(x=YourNamePCA$YourGroup)
# Kernel density estimation for group classification
Hscv1.Bd <- Hkda(x = YourNamePCA$x[,1:3], x.group = as.factor(YourNewfileName.dt$spec$YourGroup), bw = "scv", pre = "sphere")
# Perform kernel discriminant analysis
Tt.kda3.Bd <- kda(x = YourNamePCA$x[,1:3], x.group = as.factor(YourNewfileName.dt$spec$YourGroup), Hs = Hscv1.Bd)
# Plot 3D plot
plot(Tt.kda3.Bd, colors = c("YourColour","YourColour","YourColour", "etc"), drawpoints=TRUE, col.pt=c("YourColour","YourColour","YourColour", "etc"), box=FALSE, display="rgl")
# Add text
text3d(YourNamePCA$x[,1:3], texts = YourNewfileName.dt$spec$LabID, pos = 1, cex = 0.8)
######################################################################
# Plot lollipop graphs for PCA components #
######################################################################
PC1 <- plotRefToTarget(
YourNamePCA$shapes$shapes.comp1$min,
YourNamePCA$shapes$shapes.comp1$max,
method = "vector",
mag = 1,
gridPars = gridPar(pt.size = 0.25, pt.bg = "red")
)
PC2 <- plotRefToTarget(
YourNamePCA$shapes$shapes.comp2$min,
YourNamePCA$shapes$shapes.comp2$max,
method = "vector",
mag = 1,
gridPars = gridPar(pt.size = 0.25, pt.bg = "blue")
)
PC3 <- plotRefToTarget(
YourNamePCA$shapes$shapes.comp3$min,
YourNamePCA$shapes$shapes.comp3$max,
method = "vector",
mag = 1,
gridPars = gridPar(pt.size = 0.25, pt.bg = "green")
)
######################################################################
# Export PCA scores for external analysis #
######################################################################
# Export PCA scores for external analysis
write.table(YourNamePCA$x, file = "FileName4.tsv", sep = "\t")
# Set working directory
setwd("INSERT_PATH_HERE")
# Load necessary libraries
library(geomorph)
# Load data
load("filenameSlidLandmarks.bin")
load("YourNewfilename.dt.bin")
###########################################
# HYPOTHESIS 1: ECOLOGICAL ALLOMETRY TEST #
###########################################
## Create Geomorph Data Frame (GDF) ##
gdf <- geomorph.data.frame(
shape = filenameSlidLandmarks$coords,
YourGroup = filename.dt$spec$YourGroup,
ind = filename.dt$spec$LabID,
Csize = filenameSlidLandmarks$Csize
)
# Reduced (Null Hypothesis)
YourName_Reduced <- procD.lm(shape ~ log(Csize) + YourGroup, data = gdf, SS.type = "I", iter = 999)
anova(YourName_Reduced)
# Full Model (Alternative Hypothesis)
YourName_Full <- procD.lm(shape ~ log(Csize) * YourGroup, data = gdf, SS.type = "I", iter = 999)
anova(YourName_Full)
# Model Comparison
anova(YourName_Reduced, YourName_Full)
# Pairwise Comparisons
YourName_Pairwise <- pairwise(fit = YourName_Full, fit.null = YourName_Reduced, groups = gdf$YourGroup)
summary.pairwise(YourName_Pairwise, test.type = "dist", stat.table=TRUE)
summary.pairwise(YourName_Pairwise, test.type = "VC", stat.table=TRUE)
summary.pairwise(YourName_Pairwise, test.type = "DL", stat.table=TRUE)
# Plot Allometry with Prediction Line
# Set colors
color <- c("YourColour","YourColour","YourColour", "etc")
e<-as.factor(gdf$YourGroup)
color <- color[as.numeric(as.factor(e))]
# Plot allometry
plotAllometry(YourName_Full, size = gdf$Csize, logsz = TRUE, method = "PredLine", pch = 16, col = color, cex=1.5)
# Add legend
legend(x = "bottomright", legend = c("YourGroupName1", "YourGroupName2", "YourGroupName3", "etc"), cex=0.8, fill=c("YourColour","YourColour","YourColour", "etc"))
##############################################
# HYPOTHESIS 2: SEX AND ECOLOGICAL INTERACTION #
##############################################
## Filter Data for different groups (e.g., Remove Unknown Sex and Offshore groups) ##
Indices <- filenameSlidLandmarks.dt$spec$Sex != "Unknown" & filenameSlidLandmarks.dt$spec$YourGroup != "Offshore"
filenameSlidLandmarks.dt$gpa.sh <- filenameSlidLandmarks.dt$gpa.sh[, , indices, drop = FALSE]
filenameSlidLandmarks.dt$cs <- filenameSlidLandmarks.dt$cs[indices, drop = FALSE]
filtered_spec <- filenameSlidLandmarks.dt$spec[filenameSlidLandmarks.dt$spec$Sex != "Unknown" & filenameSlidLandmarks.dt$spec$YourGroup != "Offshore", ]
filenameSlidLandmarks.dt$spec <- filtered_spec
gdfCleaned <- geomorph.data.frame(
shape = filenameSlidLandmarks.dt$gpa.sh,
sex = filenameSlidLandmarks.dt$spec$Sex,
YourGroup = filenameSlidLandmarks.dt$spec$YourGroup,
ind = filenameSlidLandmarks.dt$spec$LabID,
Csize = filenameSlidLandmarks.dt$cs)
# Reduced Model (Common Allometry)
YourName_Reduced_Interaction <- procD.lm(shape ~ log(Csize) + YourGroup + sex, data = gdfCleaned, SS.type = "II", iter = 999)
anova(YourName_Reduced_Interaction)
# Full Model (Unique Allometry)
YourName_Full_Interaction <- procD.lm(shape ~ log(Csize) * sex * YourGroup, data = gdfCleaned, SS.type = "II", iter = 999)
anova(YourName_Full_Interaction)
# Model Comparison
anova(YourName_Reduced_Interaction, YourName_Full_Interaction)
# Pairwise Comparisons
YourName_Pairwise_Interaction <- pairwise(fit = YourName_Full_Interaction, fit.null = YourName_Reduced_Interaction, groups = interaction(gdfCleaned$YourGroup, gdfCleaned$sex))
summary.pairwise(YourName_Pairwise_Interaction, test.type = "dist", stat.table=TRUE)
summary.pairwise(YourName_Pairwise_Interaction, test.type = "VC", stat.table=TRUE)
summary.pairwise(YourName_Pairwise_Interaction, test.type = "DL", stat.table=TRUE)
# Plot Allometry
# Set color
color2 <- c("YourColour","YourColour","YourColour", "etc")
e<-as.factor(gdfCleaned$YourGroup)
color2 <- color2[as.numeric(as.factor(e))]
Pch1 <- c(15, 17)
f <- as.factor(gdfCleaned$sex)
Pch1 <- Pch1[as.numeric(as.factor(f))]
# Plot
plotAllometry(YourName_Full_Interaction, size = gdfCleaned$Csize, logsz = TRUE, method = "PredLine", pch = Pch1, col = color2, cex=1.5)
# Add legend
legend(x = "bottomright", legend = c("YourGroupName1", "YourGroupName2", "YourGroupName3", "etc"), cex=0.8, fill=c("YourColour","YourColour","YourColour", "etc"))
###################################
# HYPOTHESIS 3: SEXUAL ALLOMETRY #
###################################
# Reduced Model
YourName_Reduced_Sex <- procD.lm(shape ~ log(Csize) + sex, data = gdfCleaned, SS.type = "I", iter = 999)
anova(YourName_Reduced_Sex)
# Full Model
YourName_Full_Sex <- procD.lm(shape ~ log(Csize) * sex, data = gdfCleaned, SS.type = "I", iter = 999)
anova(Allometry.Sex)
# Model Comparison
anova(YourName_Reduced_Sex, YourName_Full_Sex)
# Plot Allometry
color3 <- c("YourColour","YourColour")
e<-as.factor(gdfCleaned$sex)
color3 <- color3[as.numeric(as.factor(e))]
plotAllometry(YourName_Full_Sex, size = gdfCleaned$Csize, logsz = TRUE, method = "PredLine",
pch = 16, col = color3, cex=1.5)
###############################
# COMBINE ALL PLOTS IN ONE FIGURE
###############################
par(mfrow=c(1, 3)) # 1 row, 3 columns
plotAllometry(YourName_Full, size = gdf$Csize, logsz = TRUE, method = "PredLine", pch = 16, col = color, cex=3)
plotAllometry(YourName_Full_Interaction, size = gdfCleaned$Csize, logsz = TRUE, method = "PredLine", pch = 16, col = color2, cex=3)
plotAllometry(YourName_Full_Sex, size = gdfCleaned$Csize, logsz = TRUE, method = "PredLine", pch = 16, col = color3, cex=3)
############ 3D Clouds, with Offshore ##########
# Load required libraries
library(readr)
library(ks)
library(rgl)
YourName_Full <- procD.lm(shape ~ log(Csize)*YourGroup,
data = gdf, SS.type = "I", print.progress = FALSE, iter = 999)
color <- c("YourColour","YourColour","YourColour", "etc")
e<-as.factor(gdf$YourGroup)
color <- color[as.numeric(as.factor(e))]
PlotName <- plotAllometry(YourName_Full, size = gdf$Csize, logsz = TRUE, method = "size.shape", pch = 16, col = color, cex=1.5)
legend(x = "bottomright", legend = c("YourGroupName1", "YourGroupName2", "YourGroupName3", "etc"), cex=0.8, fill=c("YourColour","YourColour","YourColour", "etc"))
# Create density clouds for the first 3 PCs of the above PCA plot
Score <- PlotName$size.shape.PCA$x[,1:3]
Labels <- factor(x=filenameSlidLandmarks.dt$spec$YourGroup)
Hscv1.GM2 <- Hkda(x = PlotName$size.shape.PCA$x[,1:3], x.group = as.factor(filenameSlidLandmarks.dt$spec$YourGroup), bw = "scv", pre = "sphere")
Tt.kda3.GM2 <- kda(x = PlotName$size.shape.PCA$x[,1:3], x.group = as.factor(filenameSlidLandmarks.dt$spec$YourGroup), Hs = Hscv1.GM2)
plot(Tt.kda3.GM2, colors = c("YourColour","YourColour","YourColour", "etc"), drawpoints=TRUE, col.pt=c("YourColour","YourColour","YourColour", "etc"), box=FALSE, display="rgl")
text3d(PlotName$size.shape.PCA$x[,1:3], texts = filenameSlidLandmarks.dt$spec$YourGroup, pos = 1, cex = 0.5)
############ 3D Clouds, excluding a group (e.g., Plot Coastal Only and remove the group offshore) ##########
## Filter Data (e.g, Remove Offshore Group) ##
Indice <- filenameSlidLandmarks.dt$spec$YourGroup != "Offshore"
filenameSlidLandmarks.dt$gpa.sh <- filenameSlidLandmarks.dt$gpa.sh[, , Indice, drop = FALSE]
filenameSlidLandmarks.dt$cs <- filenameSlidLandmarks.dt$cs[Indice, drop = FALSE]
filtered_spec <- filenameSlidLandmarks.dt$spec[filenameSlidLandmarks.dt$spec$YourGroup != "Offshore", ]
filenameSlidLandmarks.dt$spec <- filtered_spec
# Updated gdf
gdfNoOffshore <- geomorph.data.frame(shape = filenameSlidLandmarks.dt$gpa.sh, YourGroup = filenameSlidLandmarks.dt$spec$YourGroup, ind = filenameSlidLandmarks.dt$spec$LabID, Csize=filenameSlidLandmarks.dt$cs)
# Allometry
YourName_Full <- procD.lm(shape ~ log(Csize)*YourGroup,
data = gdfNoOffshore, SS.type = "I", print.progress = FALSE, iter = 999)
# 3D plot #
color4 <- c("YourColour","YourColour","YourColour", "etc")
ee<-as.factor(gdfNoOffshore$YourGroup)
color4 <- color4[as.numeric(as.factor(ee))]
PlotName <- plotAllometry(YourName_Full, size = gdfNoOffshore$Csize, logsz = TRUE, method = "size.shape", pch = 16, col = color4, cex=1.5)
legend(x = "bottomright", legend = c("YourGroupName1", "YourGroupName2", "YourGroupName3", "etc"), cex=0.8, fill=c("YourColour","YourColour","YourColour", "etc"))
# 3D CLoud plot
Score <- PlotName$size.shape.PCA$x[,1:3]
Labels <- factor(x=dt.pca.bending$YourGroup)
Hscv1.GM <- Hkda(x = PlotName$size.shape.PCA$x[,1:3], x.group = as.factor(filtered_spec$YourGroup), bw = "scv", pre = "sphere")
Tt.kda3.GM <- kda(x = PlotName$size.shape.PCA$x[,1:3], x.group = as.factor(filtered_spec$YourGroup), Hs = Hscv1.GM)
plot(Tt.kda3.GM, colors = c("YourColour","YourColour","YourColour", "etc"), drawpoints=TRUE, col.pt=c("YourColour","YourColour","YourColour", "etc"), box=FALSE, display="rgl")
text3d(PlotName$size.shape.PCA$x[,1:3], texts = filtered_spec$YourGroup, pos = 1, cex = 0.5)